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Abstract 

A new model ecosystem consisting of many interacting species is intro- 
duced. The species are connected through a random matrix with a given 
connectivity C^. It is shown that the system is organized close to a boundary 
of marginal stability in such a way that fluctuations follow power law distri- 
butions both in species abundance and their lifetimes for some slow-driving 
(immigration) regime. The connectivity and the number of species are linked 
through a scaling relation which is the one observed in real ecosystems. These 
results suggest that the basic macroscopic features of real, species-rich ecolo- 
gies might be linked with a critical state. A natural link between lognormal 
and power law distributions of species abundances is suggested. 
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The universal properties of complex biosystems have attracted the attention of 
physicists over the last decade or so. The standard approach to complex ecosystems 
is based on the classical Lotka-Volterra (LV) 5'-species model [1] : 



where {Ni}; {i — 1, S) is the population size of each species. Here ej and aij are 
constants that introduce feedback loops and interactions among different species. 
The matrix A = (cty) describes the interaction graph [1] and even a small degree of 
asymmetry leads to very complicated dynamics, although several generic properties 
have been identified. In this respect, recent studies on a related class of replicator 
equations, where an initial random graph evolves towards a highly non-random 
network, revealed remarkable self-organized features [2]. 

A classical result on randomly wired ecosystems [3] shows that a critical limit to 
the number of species S exists for a given connectivity C. Here C is the fraction of 
non-zero matrix elements. This limit is sharp (a phase transition): below the critical 
value the system is stable but it becomes unstable otherwise. The basic inequality 
is that the community will be stable if C < a/S, where cr is a given constant. Here 
the equality defines the boundaries between the stable and the unstable domains. 
Field studies show that, in general, C = aS~^'^^ with e e [0, 1/2] [4]. 

More recently some authors have suggested that complex biosystems could be 
self-organized into a critical state [5,6]. Most of these approaches based on self- 
organized criticality (SOC), such as the Bak-Sneppcn model [5], involve a time scale 
that is assumed to be very large. But the species properties of a given ecosystem 
do not change appreciably on a smaller time scale [7] and several field observations 
suggest that real ecosystems display some features characteristic of SOC states. 
In particular: (i) the analysis of time series [8] shows that the largest Lyapunov 
exponent is typically close to zero; (ii) studies on colonization in islands show that 
after a critical number of species is reached, extinction events are triggered, following 
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a power law distribution [9]; (iii) the lifetime distribution of species (defined in terms 
of local extinction in given areas) is also a power law N{T) with 9 = 1.12 [9] 

to ^ = 1.6 [10] and (iv) in species-rich ecosystems, the number of species Ns{n) with 
n individuals is also a power law Ns{n) ~ W"^, where 7 is close to one, although 
log- normal distributions with long tails are also well known [11]. But none of these 
observations (together with the C — S law) have been explained within a single 
theoretical framework. 
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Figure 1: An example of the time evolution of a single species (inset) from a pop- 
ulation of = 5000 individuals and S = 300 species in the pool. Here C = 0.5 
and 11 = 5 X 10~^. We can see a wide spectrum of fluctuations in n{t). Main flgure: 
the number of species present in the simulated ecosystem increases up to an average 
value < S >~ 125 imposed by the matrix connectivity. 



In this letter we present a new model which can be used as a simple alternative 
to the LV formulation and follows the simple (but somewhat different) approach of 
other studies on random graphs [2]. The model involves a set ©5 of S (possible) 
species and a population size A^. Individuals of different species interact through 
a matrix fl, with a pre-defined connectivity (7,^ G [0, 1]. The Qij are randomly 
chosen from a uniform distribution p{x). Each time step two individuals are taken 
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at random belonging, say, to species i and j respectively. Then if Ojj > fi^j the 
individual belonging to j is replaced, with probability 1 — by a new individual 
of species i. If Qjj = Qji — then nothing happens. This rule is repeated N 
times and these N updates define our time step. Additionally, with probability /i, 
any individual can be replaced by another individual of any of the species from 
Qg- These rules makes our model close to equation (1) when > (^'-species 
competition model [1]). So two basic reactions are allowed: (i) interaction through 
fl and (ii) random replacement by a new species from the pool. The second rule 
introduces immigration (at a rate /i): our driving force [12,13]. 

An example of the temporal dynamics of this model is shown in figure 1 , where we 
can see a wide spectrum of fiuctuation sizes, not very different from those observed 
in some natural communities [1,8]. Starting from any given initial number of species, 
this model self-organizes towards a subset of species (with average < 5" >). The 
specific subset itself changes in time. Let C* be the probabihty of having two species 
connected. Then C* and < S > are linked by a power law < S >— k{C*)''^'^'^^'^\ 
Once this < S > is reached, complex fiuctuations can arise and under some condi- 
tions (see below) power laws are observable. The analysis of the population fiuc- 
tuations typically shows, for small driving /i and large enough C-^-, a power law 
P{n) ^ n^^f^{n). Here 7=1 when ^ S* ^ 1 and ffi{n) an exponential cut- 
off. The lifetime distribution of species is also a power law: N{T) ^ T^^ and for 
N ^ S, 6 & (1,3/2). The 9 = 3/2 limit corresponds to low-level connectivities (i.e. 
random-walk-like behavior) and ^ = 1 to the situation were interactions dominate 
over immigration. It is remarkable that these limits in the lifetimes correspond to 
those reported from field studies [9,10]. 

The limiting case for // — > shows the exact hyperbolic relation and can be 
derived from a mean-field approach. It is not difficult to see that the number of 
species will change in time following: 

^^^iN-C*S, (2) 
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Figure 2: Stationary probability distribution Ps{n) obtained from (4)- (6), here n — 
10~^. For different values of S/N and C we get different stationary distributions. 
For S = 500, A'" = 50000 and different connectivities, we get Gaussian (C = 0.001), 
lognormal (C = 0.01) or power laws. A power law with a perfect scaling is shown 
for fi = 10^® and C = 0.1 (here P{n) ^ n~^). We also show a power law with 
exponential bending (black squares). 

where the first term on the right-hand side stands for the continuous introduction of 
new species at a rate fi and the second for the instability derived from interactions. 
The previous equation gives an exponential approach to criticality S{t) = \jiN — 
{piN - C*S{0))e-^**]/C* which leads at the steady state to 

S- = ^<x-l. (3) 

A mean field theory of species abundance can be derived. Let us start with the 
master equation for this model. If P(n, t) is the probability (for any species) of 
having n individuals at time t, the one-step process is described by [14]: 

dPin t) 

— ^ = rn+iP{n + l,t)+ gn-iP{n - 1, t) - (r^ + gn)P{n, t) , (4) 

5 



where the one-step transition rates are: 

^ W{n -l\n)= C\l - /^)^^ + ^{S - 1)^ , (5) 

^„ ^ W{n + l\n) = C\l- H)^^^ + f ^1 - ^) • 
Prom (5) and (6) we see that ro — and Qn — and we define g-i — and 
rjv+i = 0. The stationary distribution Ps{n) is obtained from dP{n)/dt — using 
standard methods [14] namely by writing 

Ps{n) = 9n-l9n-2...90 p^^^^ ^ 
VnTn-l-.-ri 

and using the normalization condition J2n=o Psin) — -Ps(O) + J2n>o Psin) — 1, to 
show that 

^"^'^'^^ -to[nr^ r(A*) r(n + l-.r 
Here A* = /i*(iV - 1) and u* = N + i2*{N - 1){S - 1), where jj* = /i/[(l - i^)SC*]. 
This sum takes the form of a Jacobi polynomial P^'^\x) [15], with a — —v*, 
P — X* + u* — {N + 1) and x — —1, which can itself be expressed in terms of 
gamma functions for this value of x. Therefore using q) — r{p)r{q)/r{p + g), 
the stationary, normalized solution can be written in the form [16] 

As the connectivity (or the immigration rate) increases, the distributions changes 
from Gaussian to lognormal and to power laws (fig. 2). This is expected since very 
small connectivities will make interactions irrelevant, with fiuctuations dominated 
by random events. We can simphfy (9) in the region of interest: N and S large, 
and 3> rimax where n — 1, ...,nmax- If in addition A <^ 1, which corresponds to /i 
small: /i <^ SC*/N, we get a scaling relation: 

Ps(n) = }Cn~^ ex.p{-nii*S) , (10) 
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where K, = {fi*S)^* /r{X*) in agreement with simulations and field data. These 
results suggest a simple link between a range of statistical distributions from those 
generated by multiplicative events to those related with SOC dynamics. 
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Figure 3: Fraction of remaining species at steady state as a function of S/N and 
connectivity {N = 50000, fi = 0.001). In this parameter range, hyperbolic-like 
behavior can be observed for high values of S/N. 



From this solution we can compute the stationary number of species < S >= 
{1 — P{0))S and find the relation between this quantity and the probability of having 
two species connected, C*. Under the same conditions that led to (10), we find that 
Pg(0) = {ijL*S)^* [17] which in this parameter range gives < S >= S /i* N \n{l / {/j,* S) . 
Substituting in the expression for /i* gives 



<S>-- 



ln(i^)+lnC* 



(12) 



(1 - /i)C* 

where e is given by e^-*^ = ln((l — lJ.)/fJ.) ~ hi^^^ and ^ is a constant equal to 
A'"e~^ exp(— e~^). The additional conditions under which (12) hold are A* <C e and 
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|ln(C*)| -C |ln(e)|. The latter means that C* may not be taken too small, for 
example, if ~ 10~^ then C* > 0.1 . With these values the former condition is 
satisfied if N/S is not too large. Thus (12) holds for a wide range of parameter 
values. The relation between < S > and C* is shown in figure 3 for different S/N 
values and a quite large immigration rate {/i — 10~^). In this case, the hyperbolic 
relation is only reached for large S/N and for smaller S/N values, slower decays (i. 
e. < S >^ C-^+^) are obtained [18]. 

In summary, a simple S'-species network model that shows both lognormal (even 
Gaussian) distributions as well as SOC dynamics has been introduced. Starting 
from random initial conditions, the system evolves towards a state characterized by 
some well-defined scaling laws and statistical features. High immigration, small S/N 
relations or low connectivities lead to lognormal distributions. These are replaced 
by well-defined power laws as the interactions become more relevant. In the last 
multiplicity of metastable states, the presence of a threshold in the number 
of species Sc —< S >, the scale separation between the driving when /i is small and 
the system response fully describe our system as SOC. As a consequence of critical 
dynamics several well-known field observations are recovered. The interpretation 
of this model provides a new framework for understanding how complex dynamics 
emerge in multispecies ecosystems and suggest that some well-known time series 
from natural communities would be the result of criticality instead of deterministic 
chaotic behavior. 
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